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ADS'^ract 


Tuo-dimensional istea'ly supersonic flow over an airfoil is 
considered. A hiqhly accurate approximation and a new, 
rapidly ccnverqent nuraerical procedure partly based on it 
are developed. Examples for a symmetric airfoil over a 
ranqe of Kach numbers are given. Several interesting 
features are found in tli? calculation of the tail shock and 
the flow behind the airfoil. 
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1 , In tr od u c-*-i on 

This oaoer contains an investigation of two-diroensional 
supersonic gas dynamic flows. Although the final step in 
our investigation is numerical, we use methods which 
incorporate our analytical and physical know edge of such 
flows. The approach is woU-suited both for numerical 
integration and for the interpretation of the result.ing flow 
phenon-ena. In the present investigation, several new or 
little-known effects conesrnirg the tail shock and flow 
behind a two-dimensional airfoil emerge. Discussion of 
these is reserved for section 6. A preliminary version of 
this approach for the case of one-dimensional unsteady flow 
has already bean reported (Sirovich f. Chong 1980, Chong 8 
Sirovich 1980), 

For problems in which the shock waves are weak, 
variations in entropy and one Riemann invariant are 
third-order effects, and the solution is approximately given 
by a simple wave, or Prandtl-deyer expansion (Friedrichs. 
1948; see Lighthill 1960 for corrections and extensions). 

In fact the numerical change in the F.ioraann invariant is 
significantly smaller than tiiat in entropy (section 3). 

This suggests that a larger class of flows can be viewed as 
the interaction of a simple wave and an entropy variation, 
which in turn suggests that streamlines and principal 
characteristics be used as coordinates, Adamson (1968) has 
used a similar coordinate system in another context. 


Th<^ useftilnass ol: -.his zransf oirinat ion is also related *.o 
shock expansion theory, a method £or calculatinq surface 
pressures whicn qoes back to Epstein (1931). It depends 
upon the tact that reflections of the outgoing waves on the 
principal characteristics are weak and tend to cancel each 
other (Hay«s f- Probstein I'-jriG, Eahony 1955). Shock 
oxnansion tbeo*.y has been extended to include computarion of 
the full flow field (Eaqers, Syvertson fV Kraus 1953, Eeyer 
1957), Jones (1963), in another approach, bridges shock 
expansion ‘•heory and simnlG wave theory by considering 
slowly varying perturb/ations of the latter. Our approximate 
solution is closely related to shock expansion theory. 

i'fe use ’■.he approximate solution as the first step of an 
iterative numerical method to compute the exact solution. 
This prooedure is quite distinct from current numerical 
methods for •’his type of problem. For the most part such 
approaches apply a variety of differencing schemes to the 
gas dynamic equations in their standard form (for a 
comparison of several methods sea Taylor, Kdefo 6 Masson 
1972). ImnlGmentation is then relatively simple, hut may 
require many mesh points and/or be subject to restrictive 
stability criteria. Methods which do not explicitly fit 
shock waves also tend to have difficulty with them, 
producina oscillations near or diffusing the discontinuity. 

A more Dowerfui method which fits the shock wave explicitly 
is the method (Babenko, et al, 1966, Holt 1 979), This 

is a rather unwieldy method, which is actually intended for 
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more coinplicateu probl'Tias than that dealt with here. The 
method of characteristics (see Liepinann K Poshko 1957, ch, 
12) , to which our method is more closely related, has also 
been found in practice t.o be unwieldy and time-consuminq. 

Comparisons, of computation time are difficult to make, 
because of the many variables involved, but the procedure we 
present sliould compare favorably with others available. The 
approximate sciution is probably sufficiently accurate in 
many instances, and even in the worst cases we have 
calculated only about three or four iterations are required 
to achieve an accuracy of one percent throuqhout the flow 
field. 


2 . F ormilat io n of nro bl em 

We consider the situation shown in figure 1, in which a 
uniform flow of Mach number M^,>1 is incident upon a 
wo-d iiiiensional symmetric airfoil. It is assumed that there 
are attached shocks at the leading and trailing edges, and 
that the flow remains supersonic everywhere. discuss the 

flow in the upper half plana ahead of the tail shock. The 
tail shock and the flow behind it arc treated in appendix B. 
If the airfoil ■is not syniraetric, the flow fields above and 
below it can each be computed by the methods described here, 
independently, up to the appearance of the tail shocks. 

The coordinates x and y are scaled by the airfoil length; 
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Ri' in.iiT of the 
ORiG. :ae is poor 

th'- pressure p and the density p by their upstream values Pp 
and pp ; the velocity (u,v) = (q cos B , q sin 6) and the 
sound speed a oy the uustream sound sp«6d ; and the 
entropy s, which is set to zero upstream, by the qas 
constant B, Ka consider a perfect qas with constant 
specific hears Cy = •Vt/'"!) t^nd Cp = / c^ , tor which the 
equation of state is n = p’^expf. {^~l)sl and the sound speed 
is qiven by = p/p« The calculations here were done for 
= 1.4, 'Modifications for the case of a qas with a general 
equation of state are outlined in Appendix A. 

The equations of in viscid two-dimensional steady flow ara 
conveniently wrirten in characteristic form with the entropy 
G, the f].ov; angle 6, and the Mach angle |i = sin“i (1/M) 

(where ’-'=q/a is the local dach number) as dependent 
variabloG. All other physical quantities can be obtained 
from ‘•■hese and Pornoulli’s equation 


^-1 

a 2 + q 2 = T + H 2 , 

2 2 


The Hduations of motion are (Meyer 1960, p. 273) 


ds = 0 on streamlines: — = tan 0, 
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where F (10 is -jiven bv 

P(lO = fx tan”MVAtan ja) - p, ^ = {/+ 1) / (y- 1 ) . 

The streamlines and the C-*- and ' Gharacteristics are shown 
in figure 1. The quantities r^ = B±P (ji) are called the 
Pieinanii invariants. Another useful form of (3) and (4) is 

1 dp dy 

de ± — sin 2[i -~ = 0 on C~: — = tan (e±ji) , (3) 

2Jf- p dx 

The appronriate boundary condition at the airfoil is 

tan e = £’(x) on y = f (x) . (6) 

In addition, the solution far away from the airfoil {y-»±oo) 
must approach the upstream conditions 6 = 0, , and 

s = 0. The -jumps in 6, p, and s across the shocks are 
governed by the Rankine-Huqonict conditions (Liepmann 8 
Roshko 1P57, p. 85) . If we denote the shock angle by 7j, and 
flow quantities on either side of the shock by subscripts 1 
and 2, these conditions can be written 

1 (^f- 1) tan2 (ape, ) - 1 

tan(‘h-e,) V*-H /-I 

(1 + M2) + (1+ Mf) tan2 ('M-e ) 

2 2 


tan (e^-e, ) 
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where the shoch streiiqth z ~ (Pj.-P|)/P| is qiven by 


2V- 

2 = — [T^siiiM're,) - 1]. (10) 


Equation (P) determines the jump in p, since can bo 
written as a function of M or |i usinq (1). 


3. A n ew coordin ate system 

'veross a shock wave the entropy ami one of the Sieraann 
invarian'^-s (r“ for fj positive, r+' for • negative) change 
only at third order in the shock '^trenath. This can be 
shown by laylor series expansions (or see Courant S 
Friedrichs 194d, Section 138), In figure 2 the jumps in s 
and r~ are plotted on a Loqar irtunic scale against the 
deflection angle for various Mach numbers M,. without 
loss of generality, we can set 6|=0 and take 6j_, and hence “i 
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posiiiv'?. Tk? shock strength is the same order as the 
deflection angle, so for small 8 j_, the curves approach 
straight lines of slono three. For larger values of 6 j_, the 
rates of increase of As and Ar~ tend' to drop off somewhat. 

For a weak shock wave As and Ar~ can be considered 
negligible compared to, say, Ar^, which is a first-order 
quantity. These values also provide estimates for the 
variation in s and r~ throughout the flow field. This is 
obvious in the case of s, since it is constant on 
streamlines. In fact the same is very nearly true of r~ for 
a wide range of conditions. This will be seen in connection 
with, shock expansion thaorv l.iter on. Therefore, for a weak 
shock wave s and r~ are nearly cons'^ant everywhere, and the 
solution is approximately given by , a simple wave on the C+, 
or principal, characteristics (sea section, 4) . 

A second feature of intere.st in figure 2 is that for any 
given Mach number and deflection angle the jump in r- is 
significantly smaller than that in s. At = 5 and 
82 ^ = I* 2, for example, As = 0.19 while Ar~ is only 0.03. 
Therefore, for v/eak to moderate strength shoc.k waves, the 
flow in ■^he upper half plane can be regarded as primarily an 
interaction between the simple wave and an entropy 
variation, with r- playing only a small role, 

With ^-his in mind we introduce a coordinate system (o{,p) 
consisting of the the streamlines, o( = constant, and the 
principal (C + ) characteristics, constant. By 

definition, oi and p must' satisfy 


- 8 - 
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dt^ +<X^tan e = .), jS» ^ = '■>• 


Sinc'3 any function of «l will satisfy the first equation and 
any function of ^ will satisfy tiie second/ two conditions 
will be requirad later to fix the transformations. 

In the new coordinate system (2) becoiaes 


Sp = 0, 


( 12 ) 


or 


(oV) . Similarly (3) and (4) becoae 


G3 + P(|i))^ = — sin 2ji s’ {<<) / 

2/ 


Is > ‘ , 


i.iij i'- nns 

(13) 


(1 * » *^) (e-?(|i)) 


- •— sin 2u s' (c<) , 

2# 


(14) 


w h ci r e 


= tan(9-p) 

0 Cj 4 + o(^ tan (6-|i) 


1 - tan e tan ji x^ 


(15) 


The second expression for w follows from (11) and the 
functional relation 


°<y\ 

1 

/ 

-Xp \ 

1 F* 1 



^ / 


(Id) 


Using (15) and (13), equation (14) can bo simplfied to 


{e-P(ji))^ = (1 - tan e tan |^) — . 


xj, 


(17) 


Sinct X ani y are now depcmlont variables, ^wo additional 
equations arp required for them. These are furnishel by 
(11), which by (1G) can be written as 

Vp = xp. ta'^ (e+ji) . (IB) 

Equations (12) r (13), (17), and (IB) are five equations in 
five unknowns: e, ji, s, x, and y. 

It is rossiole to eliminate y from the equations 
immcaiatelv bv settinq y^^ = y^^^ in (13). This -gives 

0 = [ tme-tan (6+ji) ] — e^sc-c^e— ■ - (e+|i) ^sec^ (e+)i) , 

which, using (17), can be written 

0 = + (j) + p (ji) )p cot (1 + (6 + |i)p tan (e + )]) . (19) 

tJio first and third terms are the ^ derivatives of log 
and -Log cos (e + |i) , respecti voly. Tho second term can also 
bo in*-pqratPd explicitly; the result is A log a, where we 
recall A= (/■+ 1 ) / (sf ~ 1 ) • Therefore (19) has as a first 
int eqral 

= A (oO a-^co3 (etji) , (2.1) 

where A (o|) is an arbitrary function to be determined later. 
Hence x can be written 

0 ( 

x(of,^) = x(0,p) t i" A (e|) a-'^cc.-3 (e + p) do{. 


( 21 ) 
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Sirailqrlv, troai (1*3) wo qo’; 

«( 

= vV'>tp + J \ (<0 a-^siri (e+)i) dflf. (22) 

O 

Sc far Of and ^ have not boon .5])''cifiod beyond saying <( = 
coiiGtaiib on streamlines and p, = constant on C+ 
characteristics, Thi boundary and shock conditior.s in the 
fl(p-Dlano can be sinplifisd somewhat by normalizing o( and p 
appeonr iatoly. We let the airfoil surface bo o( = 0, and 
nermalizo by setting p = x at o| = 0, Ihe boundary 
condi*-ion (f) then becomes 

jc(0,p) = j3 , y(0,p) = f(p), e(0,p.) = tan-if'(^). ( 23 ) 

Cno convenient way of norm.ili'/.ing d is to take the front 
shDCK ranqle ')j(o() to be given !)v 

tan '>|(oO = (l-ei)-an 'Yj{0) + d tan ji^, , (24) 

wliors ^^{O) is Known from solving the shock conditions at ':iio 
leading edge, and jij, is the uDstroam Mach angle, which the 
shock apnroachis far away from tiva airfoil, (We assume that 
'>) is a strictly decreasing function,) The flow field in the 
upper half piano thus is mapned i.nto a finite region in the 
dp- plane, as shown in figures 4 and 9. The principal 
characteristics become vertical linos, and the streamlines 
become licrizontal lines. The front shock maps into some 
curve p (of) , and the rear shock into two separate curves 
p^(c() and p 3 («f) . Flow variables on p^ (d) denote values just 
to the left of the rear shock, and thoa^a on denote -ha 


values iust to the right. The dicjcuj-’-Miou of these is la^t 
re appendix P. 

With th« shocK anql* a given Cunct.ion^ the shock 

conditions (7)- (<3) can he iiTniediato Iv solved for 
CV, ^(ci ) ) , and s {ci) . The shock ^{o() i’-solf will in general 
depend on the rest of the solution, however. 

At Bid) ^hc condition 


tan 7j 


'll - ^ 

dx (cl) 


RilPiWc'Uv ^ 
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"it 01 '' 


(25y 


must he satisfied. Usino (Id), this can be writtc-n 


+ b id) xp, = '1 on |3 = p(c^) , 


( 26 ) 


w l! e r « 


h(cO 


'“an >) - tan H 

BMoG 

' tan Tj - tan (ft + ji) 




Si)h3titu*-.ion of (21) for x In (26) aives n linear intcvjrnl 
eq nation for A (o( ) : 


«( 

A (o<) 0 (c(, p(c() ) e b((X)[1 + 5 A (CT) 0 - (cr,p(o<) ) uer] = 0 (27) 

where Q ~ a-^cos (S + p) . If the solution for 6, ji, and s is 
known in the (V^-plane, this aquation can i/o solved for A (^) , 
and tlie transformation back to the rdiysical plane coinputad 
with (21) and (22). In qeneral however the solution in the 

q(0- plane depends on x throuqh (17) . 

# 

Up to this point the equations in «<'6-cooruinates have 


been ticrivpd without: iDpcoximation, and hence are equivalent 
to the oriqinal set ('?)• (4), 


4. Aiinrox iira te solution s 

Simple wave ao p cox ima ti on 

As irie:rticned earlier, in a problem wirh weak shock waves 
dovia*-.ions In s and r~ from ♦•heir upK-^-reata values are 
third-order quantities and can be r.eqlectc-d; that is, ir can 
;icj=urnf=d thar s = *) and r- = -P(|io) everywhere. The 
solution of (2) -(4) then is a simnle v/ave, in which all 
quantities are constant on the principal characteristics, 
v.hich in turn ars straiqht lines: 

e = tan-if'(^), [1 = P-MS-t-P (po) ) / s = d 

on C-*-: y = f (p) + (x-p) tan (e + |3) . 

This approxiinaeion is due eo Friadrichs (1 948). (Friedrichs 
furthc-r simpliried the problem by noglecting terms of third 
ordar and tiqhsr throuqhout the calculation,-.) The solution 
satisfies the boundary condition, but can sa-cisfy only one 
or the three conditions at the shock. However, since two 
quantities (s and r~) are conserved up to third order across 
the shock, if one condition is satisfied the other two will 
be satisfied up to third order. It is convanient to retain 
the shock condition on 6, equation (7) , which can be solved 
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£or tht^' shock 4nql« as a function of © (^) . Denotiny tha 
shock bv [X Ip) ,Y|p) ] we have then 

KI^PiJlH'Uv'rmiTY OF THE 

Y Ip) = f ip) + [ X|p) -p]tan (©♦•ji) I A(.,) 13 POOR 

Y* |p)/X’ ip) = tan '>1 . 

Elinina^in i Y |p) between these two produces a linear 
first-order ordinary dif ferenrlal equation for X (p) , for 
which an explicit solution is easily found. 

Shock ex pansio n th eotl^ 

licca'ise simple wave theory takes s and r~ constant at 
tlioir upstream values, it can be expected to be least 
accurate near the airfoil, where the shock is stronqest and 
the deviation from upstream conditions is the greatest. An 
improved approximation in this region can be obtained using 
shock expansion theory, in which s and r~ are assumed to be 
everywhere equal to their values behind the shock at the 
leading edge, say s = Sq and r~ = r;^. This leads to a 
sli-gntly modified version of the simple wave solution, 
namelv 


e = tan-if ip) , ji = f s = Sg 

on C+: y = f|p) + |x-p) tan (e + ji) , 

This approximation produces a very accurate solution at the 
airfoil, even for flows with strong shocks, in which s and 
r- are no-^ at all constant globally. Errors should be 
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expected becauoC th^-- C~ (compression) waves which ars 
produced by reflection of the C+- (expansion) waves at th = 
front shock have b«en neglocted. Hayes & Probstein (1966) 
explain that taese reflections however are fairly weak, an-i, 
more imnor'antlv, arc nearlv cancelled by expansion wavas 
produced by tha interaction of the C+ waves with the entropy 
or vorticity layers. 'lahcny (1955) aives a similar 
explanation. 

The shock expansion solution rapidly loses accuracy as 
the distance from the airfoil, increases. This is in 
contrast to simple wave theory, which is accurate at 
infinity. 

Pr ese nt a por ox iinat io n 

We make the assumption that the flow angle 6 is 
approximately constant on principal characteristics. This 
holds true in the simple wave solution, and in general is 
closely related to s'nock expansion theory. This 
relationship will be brought out la'^er. If = 0, then 
(17) reduces to 

!«-P(li))^ = 0 or e-P(|0 = -Pot^^). (2d) 

where P^ (ci) = ?C |-i ( ^ (p/ ) ) ] - (®<, jj(e^) ) is given explicitly 

by the shock conditions. Substitution of 6 = F(n)-Po(<<) in 
the remainina aquation, (13), then gives 
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2t 


sin 2|i s'(<i<) 


(29) 


Pq (of) aiul s (of) are known functions, so (29) can be reqarded 
ns an orviinary difffcrontial oquation for p, in which p 
enters only as a parainotnr, through the initial value 


p (9,p) = + tan-if • (^) ]. 

Fqiiation (to) nonlinear, but can be rstadiiy solved using 
standard nuinerxcal methods. The solution in the e/^-plane is 
then ccinnlGted by computing e(oi,^) = P (p (of , ^) ) -P ^ (of) and 
finding ♦■h^ shock ^(d) from the computed solution |i(«f,p) and 
the values p (o(/ p («>() ) given by the shock conditions. (Since 
in oractic'i this would involve intorpolat ion , it is more 
convenient to use p {‘S, ^(oi) ) as tlie Initial value and 
integrate (29) downwards along each characteristic p= p (c<) . 
ThuJi ^ (<V) can be computed from 


f'(p(cO) = tan fP[ji (0, p(cO ) ] - Po(0)} 


by lnver*-inq f.) The solution for e, p, and s in the 

plane is independent of x and y, because (17) , the only 
equation in which x or y appears, is neglected. Tlie 
transformation back to the xy-plaiic is found by solving (27) 
for A (o() (also an easy numerical calculation) and evaluating 
the integrals (21) and (22). The solution obtained from 
this approximation will satisfy the boundary condition and 
all three shock conditions, but will satisfy (17) only 
approximately. 


This approximation is related to shock expansion theory 
in the following way. Shock expansion theory shows that r~ 
is approximatoly constant along the airfoil, and, as has 
been pointed out by ’•lahcny and Skeat (1 955) and Meyer 
(V^57), since any streamline is a potential airfoil, r~ 
should he approximately constant along each streamline. 

This is just (28) . 

In the literature this assumption is employed in various 
wav-3. If r~ = -Po (<y) , then by (17) € = 6 (p) , and hence also 
P = p(p)f can be seen from (5 + ) . Taking both e = e (p) 

and p = p(^) along with r“ = “Po(o{) overdeter mines the 
problem however, since any one of 8, p, and r~ can be 
written as a function of the other two (and s) . This was 
noted by Fjgqers, Syvertson, and Kraus (1953). In their 
generalized shock expansion method (a numerical construction 
similar to the method of characteristics) , they resolve this 
by averaging results assuininc 7 r- = -P^ (o<) and 6 = 6 (p) with 
those assuming r~ = -P,, («”) P = P (f>) • While this seems 

somewhat arbitrary, it can be shown -^hat the correct result 
in fact lies between the two (sea Hayes ft Probstein 196S, p. 
995). Meyer (1957), on the other hand, implicitly drops the 
assumption p = ,p(^), and uses the solution r~ = ~P^(c^) and 
& = e(^), which satisfies (17) exactly, but does not satisfy 
(l-l) . 

It is more consistent to approach the problem in either 
of two ways: in equation (17) assume (i) the left hand side 

or (ii) the right hand side is zero. Then solve (17) along 


with the remainiiiq equation, (13) . In case (i) , the 
solution becomes 6 = 0 (^) , P = ? (p) , and s = s (e<) . The 
function 6 (p) xs dateirinined by the boundary condition, and 
P(^) must b-3 determined by the shock' conditions. It than 
happens *-hat oyer the rear half of the airfoil, p > p(1) , 
p(^) cannot be found, slnco no data is specified on the raa 
shv.K:k. This difficulty does not arise in approach (ii) , 
which is the one we adopt. This methad requires more work, 
hut has been found to be. more accurate. It also produces 
the same solution at the airfoil as shock expansion theory. 
Additional support for this choice is lent by the fact that 
the factor multiplyinq in (17) is in general quite small 
Approach (i) has however been found useful for calculating 
the flow behind the tail shock, where the otlier method is 
inappropriate (see appendix B) . 


5. 'ujmerical. method 


Ouc approxiiiiate solution does not satisfy (17), or, 
equivalently, the C~ equation (1h), In this section wg 
present a simplo iterative method -for correcting the 
solution so that it will satisfy all the equations and 
conditions. 

The approximate solution is first computed on a 
rectanqulror grid in the plane, as shown in figures and 
9 . The same grid points' are used in the iteration scheme. 
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The front shocK p(<x) is therefore kept fixed throughout the 
itera+-ions. This fixes the normalisation of , so for every 
iteration bevoml the original approximation ^ (<r) is nor 
given by (24) and must be found as part of the solution. 

This also implies thato(=1 will no longer correspond exactly 
to Y= Off. 

Given the approxiinats solution for; 6, p, s, and x in tha 
tiC^-nlane, a corrected value of r- is couiputod from the C- 
oguation, (14), starting at the shock witn the value given 
by the shock conditions and integrating along C~ 
characteristics: 

c- 

Tlis new value r~(0,^) at the airfoil determines a new value 
of r^ there, since r+ = 2e - r~ , and 6(0,^) is given by the 
boundary condition. Kith this as an initial value, a new r+ 
can be computed everywhere by integrating (13) along C+ 
characteristics : 

/ I 

= r*(f),p) t j ^1’ l«) d"*- (31) 

o 

The solution given by r~, r+, and a will satisfy (13) , (14) , 
and the boundary condition. However, the new value of 

p (c/) ) will no^ in general satisfy the shock conditions, 
and hence will imply new values of 0j(e<) , s (o() , and 
( 4, p(p?) ) . This furnishes a new initial value for 
in*-eqrating (3.)) , which is used to start the next iteration. 
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In siirniiiary, the steps of the method are as follows: 

I. Compute approximation^ fixinq front shock and grid 


in plane. 

II. Compute r-(fl|,^{«<)) from Iht shock conditions and 
evervwhcr'’ from (3.)) . 

III. Compute new r+(0,p) = 2« (9,^) -r- (0 ,^) with 6(0,^) 
from the boundarv con«li*:io!' ind r~('j,^) from step 
II. Coraoute r+(e<,|3) ''vor vvluir^' from (31). 

IV. Compute new values B = (r++r-)/2 and |i = 

P“i { (r+-r~) /3) G vorywliCro. 'v'ith the new solution 
at the shock compute a now 9j(c<) and s (o<) from the 
shock condi+- io ns. ^Iso recompute A (d) and x {d r ^) • 

V. Checi\- for convertrence and either gc to step II or 
stop and compute y(<^,^). 

Tills scheme has been implemented using second-order 
numerical methods (trapezoidal rule, etc.). The details are 
given in appendix C. 

The evalucuticn of the integral (30) at step II is the 
most complicated calculation in the' p roced ure. In place of 
eguvUtion (14) we might have integrated (17), which is 
equival^n'- and has the advantage that r~ is differentiated 
only with respect to p. In practice however, this does not 
work very well. The itoration procedure does not converge 
as quickly, and may not converge at all without modification 
(see Chong girovich 1980). Ke attribute this to the fact 
that small variations in r~ are naturally propagated along 


tlio C~ characteristics. 


(•, nn 1 nif^cussion 


.i j ! 


-fo 


The oxecution of tho nnaiericil procedures which we have 
discussed is both rapid and inexpensive. As a rcsul-. we 
have perforned a nuraber on calculations for several airfoils 
over a ranqe of flach nuiobe-rs. The results presented in the 
fiqures are for a synun^tric circular arc airfoil with 
thickness ratio 0,?5, at upstream dach numbers li,, = 2.5 and 
7.5, In seme fiqures the intermediate case Mq = 4,0 is also 
shown. These cases were chosen in part for the intersstinq 
effects they exhibit. 

The i-^eration scheme is foinnl to cenverqo quite rapidlv, 
based on a comparison of th'= solutions at succcsive 
iterations. In the followinq table, the maxima (over all 
qria roiiits) or the diffarenc^^s in tlue values of e, ji, and x 
are qiven for the case Mo=7.5, 



xteration 


1 



A 0 /e (d,0) 

Ap/jl 

A x/x 

0.0595 

0.1170 

0.3701 

1) . 0 1 4 1 

0,n096 

0.1221 

0 . )ooa 

0.0007 

0.0112 

0.0002 

0.0006 

0.0017 

0. 0001 

0.0003 

0.0009 


The greatest differences are in x and occur within a few 
qrld points of o{ = 1, where x-»oo, TliC errors in x are smalle 
closer to the airfoil. For thinner airfoils and/or lower 
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^lach nujnbfics, fewer iterations are reguireu for the same 
accuracy. In the case of a 1)^ thick parabolic arc airfoil, 
for exaniple, even at m„ = 10 the difference between the 
approwisate and exact solutions is less than one percent in 
e arm ji and six per cent in x. In such a case there is 
little reason to qo bevond the approximate solution. 

The case .Mp = 2.5 is discussed in Holt (1977). Figure 3 
contains a ccraparison of the leading shock when computed by 
our method and the BVLd method. The small difference is 
probably attributable to coDying errors. In addition wa 
indicate our approximate theorv and that which emerges from 
a variant of snock expansion theory (the latter and the BVLF 
curve are taken from Holt, p.77). In this case, our 
approximate solution is indistinguishable from the exact 
solution. 

■ie continue the presentation of this case by indicating 
in fiqure h a portion of tne arid in the t^p-plane used in 
the numerical intearation, (The increment between the 
streamlines shown is 0.1; the computations were done with 
^o( = d,d2‘^.) The triangular region corresponds to the flow 
behind ‘-be tail shock (see appendix B) , Figure 5 gives the 
entropy s (ct) in the region between the front and tail 
shocks, and s^lcH) in the region behind the tail shock. The 
deflection angle e is plotted versus d in figure 6 on each 
of the C+ characteristics shown in figure 4. A 
specification of B in ♦■he region behind tiie tail shock is 
deemed unnecessary since it is very nearly zero everywhere 
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{ i B I <0. .)05) . Fiqure 7 contains ths variation of r~ versus 
on each streainiine of figure 4, Finally, in figure- *3 the 
streamlines, characteristics, and shock wavns are shown 
transformed back to the xy-plane. Those figures completely 
determine the solution. In particular, corresponding to any 
point (x,y) of the f’hysical plane we can doterinins the 
coordinates The entropy s is then found from figure 

5 and the dcflaction angle <3 from figure tj. The latter and 
r~ from figure 7 determine the Prandtl angle P(|i). All 
o*hor flew properties then follow. 

For comparison w<-- give in figures 9-13 similar data for 
the same profile at ;i^ = 7.5. In this case |0|<O.O2 behind 
the tail shock. 

'.be Login our discussion of these results Ly considering 
the limit x to . Far behind the airfoil the pressure 
becomes constant, p -> 1 , and as a consequence 0 0. It 
then follows from the etjuation of state that 

a2 = expC- (/-1)S3 (o<) 

where s j (of) is given either by figure 5 or figure 10. From 
(1), we can then compute the velocity q at infinity. This 
is shewn in figure 14 for Ho = 2.5, 7.5, and the 
intormediare case Iq = 4.0 (same profile). As a result of 
the non-uniform entropy, the flow at infinity has a 
vorticity distribution. 

Another feature cf interest is the entropy variation 
along the tail shock (figures 5 and 10). This has a 


two-scalc anpea raticp, especially at the higher Mach number, 
which shews a very rapid decreasa in s-rnngth in the initial 
pertion of the shock. The slower variation in entropy 
follows that induced bv the front shock. Looking at figure 
13, we see that the streairlines spread apart rapidly as the 
flow passes the inidchord position. The flovf inclination at 
the tail shock thec<-£oi.e decreases rapidly, which results in 
a corresponding decrease in shock strength. 

iincther important effect is also at work in this region. 
The gas, which is ccni pressed at the front shock, in 
followinq the profile past the niidchord experiences a rapid 
expan«^ion, which is sufficiently strong so that the local 
Macn number at the trailing edge exceeds the upstraani value 
(K = S.D3 for the = 7»5 case). Sines fi is small, the 
larje negative value of e on the after part of the airfoil 
causes -he principal characteristics to have negative 
slopes, so that waves originating there must intersect the 
tail shock near the ai.rfoil. This quickly cuts off the 
recovery process. /<s a result the ilach number along the 
tail shock falit- off rapidly, which augments the rapid 
decrease in strength of the tail shock. For the case 
Hj, = 7.5 the -lach number along the shock even falls below 
7.5. 

A feature which is somewhat difficult to perceive from 
figures 8 and 13 is that the tail shock angle is not 
monotcnic. in figure 15 the variation of the slope of tho 
tail shock is given for the three cases wc have discussed. 


For n;»ch c-Tae the «^hool< angle decroascs on leavina tho 
trailim c>1qc. (This result- has been verified inclepondentlv 
by Jane? C. Townsorul 1^7b, nsin^ a different numericil 
raetnod developed bv 'lannel n. dala.-i.) Tliis is contrary to 
what is observed for lower Mach numbers or thinner bodies. 

We have seen that the angle of the incident flow decreases 
along the shocK. If the Mach number to the left of the 
shocK were constant, this would predict a decrease in shock 
•angle. The decrease in Mach number along the shock tends to 
liave the opposite otfect however, to increase the sliock 
angle. In these cases, near the trailing edge the 
dccrecasi.nq flow ancle domindtes. For fign Mach numbers the 
shock angle is more deneni'i-viit the flew angle ‘•han on the 

Mach number, ris can be seen from the fact that the shock 
polars for different Mach numbers aoproach a limiting curve 
as M eo (see j.g, Lienmanu & <"oshko 1957, p. 87). For 
lower Mach numaars or thinner airfoils the affect of 
decioasinq Maca number dominatos the effect of decreasing 
flow angle. 

heturninq to figure 15 we also see that fo;e =7.5 tho 
slcce undergoes a second oscillation in wiiich it rises above 
the Mach angle at infinity. This is explained by the rapid 
fall-off of Mach number along the shock, below its valua at 
infinity, A final item of note in figure 15 is that for 
Mo = 7.5 *he shock slope actually starts off with a value 
which is greater than at infinity. The expansion process 
along the profile nroduces a relatively high Mach number at 
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. ui 

'■\i POOH 

•thr: nrailinq eiqo, bu* it is not qc^H. t enough to bring the 
rloptt below the une-renin value, 

Tht: nreesura dieribution behind bhe airfoil is 
interestim. in figures 1f> and 17 wc show the pressure 
contours behind the tail shock and the values ot loq p on 
the x-avis for x>1 , for the cases Mp = 2,5 and 7. 5, The 
pressure on ’•ho tear part of the airfoil is so low that, in 
spite of the hiqh shock strenqth at tho trailing edge, the 
pressure 'jump t'nrouah the shock doss not bring p, up to the 
equilibrium pressure p= 1. There is a rapid pro-ssurc 
increase iir mediately behind the trailing edge, in which P 
increases above the equilibrium value, reaching a maximum 
abou*: one chord length ovit, riie return to equilibrium from 
this point is very gradual, Tlio total variation in pressuro 
behind the tail shocks is quite snail compared to that along 
*-ha airfoil surface (in terms of lou p, about 3‘. at Mo = 2.5 
and at '!<,-7. 5) . 

several other features of the calculations also merit 
mention, ■'lost notable oerhaps is ^he surprising constancy 
of the deflection angle on the principal characteristics. 

The variation in B is almost undatsctable in figure 6 and, 
as fiaure 11 indicates, the most serious departure occurs 
nc-ar the trailing edge, where it is about 10/o in the worst 
case. As mentioned already, e = 0 is an excellent 
approximation throughout the region behind the tail shock, 

A related ohencirenor is the n'sar-straightuess of tho 
principal characteristics. This however does not remain 
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nrus in -hhe rcqion behinJ *:he airfoil. 

I'Je are also in a position to examine tho basic assumption 
of shock expansion -heory, that r~ is constant on 
str samlines. i'rom equation (17) we see that this condition 
is closely rela+ ;d to the constancy of 9 on principal 
characteristics. As both figures 7 and 12 indicate, this is 
a reasonable assumption, although somew}iat remarkably it is 
better at the airfoil than in its neighborhood. In both 
cases this assumption is far less satisfaccory behind the 
tail shock. The rapid down stroke of the r~ curves also 
indicates a large value of although 6 itself remains 

quite small. 


This work was supported by the National Aeronautics and 
Space Administration under NASA Grant no. NSG 1617. 


A pp gndi x_ A: „ ca se . of. a n arbitrar y ma s 

For an arbitrary gas, the equations of motion in 
clidracter istic form can be written (Hayes & Probstein 1966, 
p , 4 3 4 ) 

ds = 0 on dy/dx = tan 6 (A1) 

de ± i dp = 0 on dy/dx = tan(et|i) (A2) 


where J = Fo / (Po ^ P) • consider i to be a 


V 


.■ 
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function of p and s. By introducing the variables 

u)(n »s) = J$(P»s)dp and A (p , s) - (p»s) /ds ^ 

which ar<? defined so that dtJ = $dp +’Ads, (A2) can be 
wrirton 


hil<] 


i t/K 


de ± dt» = ± Ads on dy/dx = tan (6±p) 


(A3) 


If 0 and A-aro now regarded as functions of p and s, (Al) 
and (a3) are tiiree equations in three unknowns; -6, p, and s. 
Equations (3) and (4) are a special case of (A3) in which 
= P (ji) and A= (sin 2p)/2/. 

The transformation to -coordinates goes through for the 
most part as before, Equarions (12)- (14) in the general 
case become 


- 0 

(6+u) )^ = As» («0 

d 0 

9o( 7)6 

where w is still given by (15). The counterpart of (19) is 

^ciB 

0 = + {p + (4))p cot p + (B + p)p tan (e + p) 

^ci 

This equation can in principlf^ be solved in the same manner 
as (1*^), but we do not have an explicit integral of the 
simple form (2d) . 

The assumption = 0 in the general case implies 
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of 

. ■ . " IF 

> 

(6- t«) ) - 0 or fe-tJ= "l»^oi‘<) • resultinq approxiniation 

r 

can fco 9xpec*-ed to bo valid at leant in cases in which the 
behavior of the qas does not differ too greatly from that of 
a perfect gas with constant specific heats and = 1,4. It 
has been shown (see Hayes t Probstein 1966, ^7.2) ’■hat shock 
expansion theory tends to lose accuracy if "f is allowed to 
approach 1, 

Appendi x B: Tail shock for a symmetr ic a irfoil 

In the general case, the solutions in the upper and lower 
Jialf planes can be computed indepsjidantly , up to the 
appearance of the tail shocks. The flows from the top and 
bottom interact behind -he airfoil, which complicates the 
computation cf the tail shocks and the flow behind them. 

The upper and lower regions behind the airfoil are separated 
by i contact discontinuity, or slipstream, whose location is 
unknown a priori. Across the slipstream e and p are 
continuous, but the other variables jump. This complication 
does not arise in the case of an airfoil symmetric with 
resDset to the x-axis; the slipstream coincides with the 
x-axis, and can be considered a rigid boundary. The problem 
is still quite different than the front shock problem, 
because the flow upstream of ■^he tail shock is not uniform. 

The transformation to -coordinates behind the tail 
shock can be chosen differently than that ahead of it. In 


particular. It is more proper to regard the C** 
characteristics as tlie principal characteristics, since the 
C+ waves arc only produced as reflections of the C- waves, 
wl'iich originate at the tail shock (see figure 1). On the 
other hand, we, have found that for the iteration procedure 
it is better to take the C+ characteristics as the. 
^-coordina'^es, because this has the effect of putting more 
points near the traiiinn edge, where the most. rapid 
variation in the solution occurs. However, the -approximate 
solution derived below is more accurate if the C~ 
characteristics are taken as p= constant. ?lodif ications 
for the use of the C- characteristics as coordinates in 
place o,f the C'*’ characteristics are straightforward. 

The normalizations of d and p behind the tail shock are 
also different than those ahead of it. It. is natural to 
keen o{ constant on streamline.s as they cross the shock. 

Also, rather than set x(0,p), we can normalize p so that 

the- region x>1 is mapped into a finite region in the 
ofp-plane. This was done by setting p 3 («=j) = 1 + °i/2, 
pro-ducinq the configurations shown in figures 4 and 9, 

The calculation of the tail shock can be done as 

follows. We assume that the solutions for e{«(,^), |i(«f,p), 
and x(«(,p) are known in the neighborhood of {<<) . Equation 
(26) holds at as well as at |3(fl(), and since in the 

former case x^ and x^ are known functions, we can write (26) 
as an ordinary differential equation for ^^(of): 
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tan “Jji - t'in(etn) 




= IMW, 6z(c() ) 


(01) 



tan - 

tan e Xp ^ = 

Pa (of) 

As 

yet 

'>|5_(o{) is undetermined. If 

we specify one 

variable , 

sa y 

B, 

gus*- +-0 the right of the s 

hock [ i.e, at p 

= pii<X) 1 

the 

n w( 

; can solve the shock conditions for 

and solve 


(Hi). This enables us to set up an iteration procedure 
similar to that used for the front flow. 

The approximate solution used £or *ihe flow over the 
airfoil cannot be conveniently employed for the flow behind 
the tail shock, because the nonuniform flow to its left 
makes it impossible to calculate and s(c<) a priori for 

use in (29) . Tkerefor‘= we use the simpler of the 
approximations given in section 4: B = 6^ (^), p = (p) , 

and s = Sj {«() . All the characteristics intersect the 
x-axis, where e = 0 , so (^) = 0 , and hence e = 0 
everywhere. In particular e = 0 a - ^ (^) . We can 

th.erefore solve for p?.(flf) as in the previous paragraph, and 
also compute ail other quantities on both sides of the 
shock. This determines Pj (^) and 

The computation of the transformation frouic*'^- to 
xy-coordina tes is also somewhat different. Denoting the 
known value x(o(,^j(o{)) by we must have 

X (<=«, (e() ) (e<) , which iiiiplios 

x’(«) = x^(W,^3(ot)) + pj' {«) Xp(«, ^3 (c<) ) . ( 132 ) 

An equation of the form ( 26 ) also holds at & = : 


■ 


u 


y 


i. 


>1^ 


• 'HIE 

yj<J& 
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0 = x^(^, ^,(«{)) + {*<) x^(«f, jSj C«0 ) 


(!33) 


If WG cllrninatG bp.twGen theso two, substitutG 

X(^ - A j (,l) Q (<t , p) [where we recall Q=a~^cos (6 + ji) ], anti solve 

for A 3 (oj) , W(? qpt 


A , W 


x^(w) tan'JjjU) - tan. 9 


Q(«</ ^j(w) ) tan ( 6 +ji) - tan 6 




(B4) 


from which we can determine 


x(rf,p) =x^(c^3(^)) + j A3 «y)Q (<y,^) dcf 

>3 (/») 

where o{^{^) denotes the inverse of p^Cw), A similar 
equation fol'lcws for y («V , p) , 


(35) 


The iteration scheme proceeds essentially as before. 

Given r - («'/ ^3 («<) ) - rom the shock conditions, ws integi:at 9 

(30) along C~ characteristics down to the slipstrea.in c^ = 0. 
Then we reset r+ (0 ,^) = -r- ( 0 ,^) , and integrate (31) upwards 
to i<=i) • the new i"*' and define a new 6 (a? , ^3 («<) ) , which 
is used to solve for a new shock (<<') and new functions 
'>|z{«<)/ r and r~ (<f, Pj (o() ) , with which we start the next 

iteration. 


App endix .. C _: . Detai ls o f ._n uni erical me tho d 

This appendix contains the details of each step of the 
iteration scheme described in section 5, as presently 
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iinpleirent-s-l. 

S tqp ^ T. w6 take a unitorai qrid in o< with N + 1 points 
Qfj = iA<< (i=;) , 1 , . . . , N) » where iiof = 1/M, Over the I'ronr half 
of rhe airfoil the ^ nrid points are taken ro be = jS(‘<j) 

( j=n, 1, . . , , N) , where ^{ct) is the equation of the shock as 
qiven bv the approximate solution. Over the rear half of 
the airfoil the jS qrifl noints are taken evenly spaced. 

The aoproximate solution in the -plane requires the 
solution of the ordinary diftepential equation (29) for each 
Pj . This was done by the improved Euler method, which for 
the qoneral equation dj/dt = is qiven by 

At 

Ikfi = Ik ^ ) ] (Cl) 

The approximate solution is completed by solving the 
integral equation (27) for A («*) and evaluating x and y from 
(21) and (22). Equation (27) is linear, and can be readily 
solved bv the trapezoidal rule; the integral is 
approximated by the appropriate sura, and the resulrinq 
equation is solved for A {<(_,• ) ia terms of A (cy„ ) , A (o(, ),,.., 

). The derivatives appearing in the 

integral are evaluated numerically using a three point 
scheme on the uneven ^ mesh. The integrals for x and y, 

(21) and (22), are computed by the trapezoidal rule. 

Ste p II. We can write (19) as 

I 

dr- = - ^ sin 2)i ds on C~: dfl/d<v = 


w 


(C2) 
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where w is qiven by (15). In computing v, must be 
evaluated numerically, but x^ is given analytically by (20), 
The C- characteristic through the point («{,|3) is drawn back 
to intersect the line segment between («<,p-Ap) and («l + ^«^,p) 

A A 

at the point = (o( + VAo( ) • To find this point 

A 

the equaticn d| 3 /dfl( = w is integrated .from c( to cf, using the 
improved Euler method (Cl) again, with the modification that 

A ^ 

a first approximation to w must, be found by 

interpolation. Re take 

= C 1 - 

w = V d^) + V w (q('+ ,^) V r 

and then redefine V by 

= [ 1 - 5 (W i«,p) +W) ^ 

This value is correct to second order and is used to 
interpolate the value 

= (1-T>) r- («{,p-6p) + V r^ (<if+ do{,p) 

and, similarly, p and s. The value of r~ is then 

computed by the tcape7;oidal rule: 

A ■' A A 

c~(«f,p) = r- - — [sin 2ji(*(,6) + sin 2ji ][ s («^ ) -s ]. 

4i^ 

A somewhat simpler .method can be used if (e<,j 3 ) is two or 
more arid points from the shock. We write (14) as 
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; ^~c( * ~ ^1* PMef) , 

i 

t 

^ evaluating w and at [ci,^) , and replacing r~ and by 

three-point one-sided differences employing points upward 
I and to the left of (ef,^), respectively. Then we solve for 

(«fi ^) • 

Step I I T. Equation (13) can be integrated directly. 
Again using the trapezoidal rule, we have 


r+ («f + A»l,j3) = r+(«(,^) + 

I 

^[sin 2p(o{,^) 


(C3) 

t sin 2p (<rf + Aol,^) ][ s (flf + 6oi) -3 (fi() ]. 


3 1 c 0 IV . While it is possible in theory to solve the 
sliock relations given any one quantity just to the right of 
the shock, in practice it is difficult given r'*- or r~. 

. Therefore we use the new 8 (a , p) computed from r+ and r~ , and 

; solve the shocK relations for , s (e^) , and \i(o(,pic()) as 

functions of 6 {<( , p (<j») ) . Then r~ (o(, ^ (<r) ) is redefined from 

[ these values no srart step II, 

I 

; The computations of ^ (ef) and x(o(,j 5 ) are done in the same 

way as in the approximate solution. It is not necessary to 
compute y(o|,p) until the final solution is obtained, since 
only X appears in the aquations. 

r 

i st ep V . To check for convergence, we compare 0, p, and x 

[ at all grid points between successive iterations, and stop 

when the maximum difference is less than some given 


tolerance 
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Figu re cap tlonr. 


tiqucf- 1. Supersonic flow past an airfoil: heavy lines 

.ieiio<-e shoch waves; solid linss ueno':e streamlines and 
C^- charvactafistics; dashed lines denote C- 
c ha r act eristics. 

Figure 2. Jumps in entropy s and Fieroann invariant r" across 
a shock wave as a function of 6^ (B, = 0) at various 
Mach nuirhers, , As; /Ar~, 

Figure 3. Front sheet for flow parSt a 2f.i circular arc 

profile at Mo = 2.5. , DVLR method; — , shock 

expansion mothod (both from Holt, 1977); - , presont 

approximate and sxact methods. 

Fiqure 4. M„=2, 5. Flow field in <vp-plnne. 

Figure 5. 'i^ = 2. 5. Fntropy s (c() ( ) in region between 

front and tail shocks, ami Sj (c<) ( ) in region behind 

tail shock. 

Figure 6. M„ = 2. 5. Flow angle 6 vs. o( on each C + 
characteriotic of figurs 4. 

Figure 7, ^^ = 2.5. Rieroann invariant r~ vs. ^ on each 
streamline of figure 4. 

Figure 3. '3^=2. 5. Flow field in xy-plane, 'Ihe streamlines 
and c+ characteristics corresnond to the lines o< = 
constant and p= constant, respectively, in figure 4. 

Figure 9. Mo=7.5. Flow field in <V^-plans. 

Figure 10. Mo= 7.5. Entropy s (oj) ) in region between 

front and tail shocks, and (c() { ) in region behind 


- 3H - 

tail shock. 

Fiqucp 11. 'lo=7.5. Flow angle 6 vs. a on each C*- 
charactsriuJtic of figuiro 

Figure 12. '^o=7.5. Fieinann invariant r~ vs. yS on each 
streamline of figure 9. 

Figure 13. t1 -=7.5. Flow field in xy-plane. Tho streamlines 
and C+ characteristics corresnond to the lines o( = 
constar.”; and ^ - constant, respectively, in figure 9, 

Figure 14. Velocity profiles far behind airfoil for = 
2.5, 4.0, and 7.5. Dashed lines denote asymptotic 
values do. 

Figure 15. Tail shock slope tan 7)1 vs. d for M© = 2.5, 4.0, 
and 7.5. Dashed lines are asymptotic values, tan ji^. 

Figure 16.<Ho=2.5. Upper; pressure contours behind tail 
shock; log p = -0.07(0.01)0 and 0(0.001)0.01. Lower: 
log p vs. X on x-axis. 

Figure 17. '1o=7.5. Upper; pressure contours behind tail 

shock; log p = -0.5 (0.1) 0.1 and 0. 1 (0.02) 0.28. Lower: 
log D vs. X on x-axis. 


























